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Abstract 

By fitting the frequencies of simultaneous lower and upper kilohertz Quasi-Periodic Oscillations 
(kHz QPOs) in two prototype neutron star QPO sources (4U 1636-536 and Sco X-l), we test the 
predictive power of all currently proposed QPO models. Models predict either a linear, power-law 
or any other relationship between the two frequencies. We found that for plausible neutron star 
parameters (mass and angular momentum), no model can satisfactorily reproduce the data, leading 
to very large chi-squared values in our fittings. Both for 4U 1636-53 and Sco X-l, this is largely due to 
the fact that the data significantly differ from a linear relationship. Some models perform relatively 
better but still have their own problems. Such a detailed comparison of data with models shall enable 
to identify routes for improving those models further. 
Subject headings: accretion, accretion discs, stars: neutron, X-rays: stars 



1. INTRODUCTION 

The launch of the X-ray timing satellite, Rossi X-ray- 
Timing Explorer (RXTE), led to the discovery of kilo- 
hertz Quasi Periodic Oscillations (kHz QPOs) of low 
mass X-ray binaries (LMXBs) in their X-ray lightcurves. 
The frequencies of kHz QPOs range from a few hun- 
dreds to about one thousand Hz; its time-scale corre- 
sponds to the dynamical time of the innermost regions of 
the accretion flow. Thus such signals may carry crucial 
information about the central neutron star (NS), such 
as the mass, spin frequency, angular momentum, radius, 
magnetic fields and so on. Usually the twin kHz QPOs 
appear simultaneously and the lower and upper QPOs 
are almost directly proportional with each other (see e.g. 
Ivan der Klisl[200l . 

Various theoretical models have been proposed to ac- 
count for the kHz QPO signals. Table Q] shows all the 
present models we collect. Although each model achieves 
its success to a certain extent, the origin of kHz QPOs 
is still highly debated. Furthermore many new mod- 
els emerge in recent years, and systematic comparisons 
among them have not been well studied. 

In view of this, we investigate systematically the pre- 
dictive ability of the present kHz QPO models. We focus 
on those that predict the frequency relation of twin kHz 
QPOs. The moving hot spots model is not included. 
This model performs a 3D MHD simulations of the ac- 
cretion around a NS. The simulations show that the mov- 
ing hot spots on the surface of a NS can develop oscil- 
lations in the lightcurves. However, this model does not 
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provide any analyti c relation between the twin QPOs 
(jBachetti et~aill2010l) . 

In this work, we measure the frequency relations of 
the twin kHz QPOs for 4U 1636-53 and Sco X-l, then 
fit the models with the measured results. We choose 
these two NS systems for several reasons. Firstly, both 
of them have strong kHz QPOs over a wide frequency 
range; both of them have been observed more than 10 
years with RXTE. Thus the bias from the sample selec- 
tion is minimized. Moreover, the different properties of 
these two sources allow us to discuss the predictive abil- 
ity of the models. They are typical Atoll and Z source, 
respective ly. The putative spin frequency for 4U 1636-53 
is 581 Hz (Strohmayer 2 0011 ; IStrohmaver and Markwardtl 
2002), whereas the spin frequency of Sco X-l remains un- 
known. 

In the following, we firstly describe the data reduction 
procedure. Then we fit the frequency relations with all 
the available models. The predictions of NS properties 
in each model will be presented. Finally we discuss and 



TABLE 1 

Present models. 



Models 


Reference 


Sonic-point and spin-resonance 


[1, 2, 3] 


Orbital resonance (3 models) 


[4, 5] 


Precession (3 models) 


[6, 7, 8] 


Deformed-disk oscillation 


[9] 


'-lr, -2v' resonance 


[10, 11] 


Higher-order nonlinearity 


[12] 


Tidal disruption 


[13] 


Rayleigh- Taylor gravity wave 


[14, 15] 


MHD Alfven wave oscillation 


[16] 


MHD 


[17] 


Moving hot spots 


[18] 



[11 IMiller et all jl99l); [21 ILamb and Millerl j200lT); [3] 
Lamb and Milled (2003); [4] Klu zniak and Abram owicz (20 p); [51 
Abramowicz et all (120031) ; [61 Stella a nd Vietril 1 119991); 171 IBursal 
(2005); [8] Stuchllk et al (2007); [9] Kato (2001); [10] Torok et al 
m07i); [111 IBakala et a l (2008); [12] Mukhopadhyay (200|); [13] 
Germana et al (2009); [14] Oshcrovic h and Titarchukl in999f> ; [15] 
Titarchukl 1120031); 11 61 IZhanel (120041 ); I17| IShi and Lil 120W ); [18] 
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conclude our investigation results. 

2. DATA ANALYSIS 

We have retrieved all the public archival data of the 
two sources with the Proportional Counter Array (PCA) 
on board RXTE. The observation time is from Feb. 28th, 
1996 to Sep. 25th, 2007 for 4U 1636-53, and from May 
5th, 1996 to Feb. 4th, 2006 for Sco X-l. 

For 4U 1636-53, we use the event mode data for 1156 
ObsIDs, with time resolution better than 256 [is and 
an energy band of 2-40 keV . With the similar analy- 
sis pro cedure in iBarret et all (|2006D and iBoutelier et all 
(2009aj), the PDS as well as the QPO parameters (peak 
frequency, width and amplitude) in each ObsID are ob- 
tained. For the ObsIDs with PDS containing the lower 
QPO, we trac k its time evolutio n in every 128 s. Fol- 
lowing that in IBarret et~al (|2005l ). all the 128 s PDS are 
aligned in every 30 Hz inter val of the lower QP O with 
the shift-and-add technique (jMendez et allfl998h . Then 
we search for twin QPOs in each interval. For the ObsID 
with PDS containing only the upper QPOs, we directly 
align the PDS of each ObsID in every 30 Hz interval of 
upper QPO. Again the twin QPOs in each interval are 
searched. Finally, the two parts of results are combined 
and we obtain the frequency relation. 

For Sco X-l, we analyze the Generic Binned mode data 
in 187 ObsIDs with time resolution better than 256 /is 
and an energy band of 2-40 keV. Considering the ef- 
fects of the deadtime, we use a model of two Lorentzians 
plus a power law to fit each PDS. The powerlaw compo- 
nent denotes the deadtime-modified Poisson noise; the 
Lorentzians account for the contribution of the twin kHz 
QPOs. We then apply the shift-and-add technique to 
the ObsID averaged PDS as described above on the up- 
per QPOs' frequency, because the span and significance 
of the upper QPOs are larger than that of the lower 
QPOs. The inter val of shift-and-add is 50 Hz. Simi- 
lar to the result in Mende z and van der Klisl (2000), our 
frequency relation shows some subtle structure when the 
lower QPO is around 800 Hz. 

3. COMPARISONS BETWEEN MODELS AND DATA 

In the following, we restrict the NS parameters M £ 
[1.4,2.4] Mq and j £ [0,0.3] (j = Jc/GM 2 ) in our 
fittings, corresponding to reasonable equations of state 
(EOS) of NSs (Lattimer & Prakash 2007), where M and 
j are the mass and dimensionless angular momentum 
parameter of a NS, respectively. The fitting results are 
summarized in Table 2. For some models, the best fitting 
values of M and j approach the upper or lower limits. In 
each of these cases, we made an extended fitting to relax 
the limits to M £ [1.0,4.0] M Q and j £ [0,0.5]; these 
results are presented in Table 3 and discussed in Section 

H 

3.1. Comparison with the sonic-point and 
spin-resonance model 

The sonic-point and spin-resonan ce model (jMiller et all 
ll998HLamb and Millerll200lLl2003f ) attributes the forma- 
tion of the twin kHz QPOs to the interaction between 
the orbital motion of the flow and the NS rotation. The 
interaction happens at the 'sonic point' where the radial 
inflow becomes supersonic. In the frame of the model, 
the X-ray source is a NS with a surface magnetic field 



about 10 7 ~ 10 10 G and a spin of a few hundreds Hz, 
which accretes gas via a Keplerian disk. At the sonic 
point r sp , some of the accreting gas is channeled by the 
magnetic field and then impacts the NS surface to pro- 
duce the lower QPO. Some remains in clumps with the 
Keplerian disk flow, producing the upper QPO. There- 
fore the upper frequency vi is the Keplerian frequency 
vk at r sp ; the lower one v\ is the beat frequency between 
z>k and the NS spin i/ s , i. e. v\ ~ vy, = vk — v & . The 
first version of the model (jMiller et all 1199 8) leads to a 
constan t peak separation Az/, close to v s . The second 
version (jLamb and Miller 2003) introduced inward drifts 
of gas to make Av dependent on v\ (or iaj). The inward 
drifts make v\ greater than vb and vi less than z>k for a 
prograde gas flow, 

v\ « vb/(1 -Vd/v B ), (1) 

V2 « v K (l - -Ucl/Ug), (2) 

where v c \ is the inward radial velocity of clumps near 
r sp , v g the characteristic inward radial velocity of gas. 
v c \ and v s are supposed to be approximately constant 
d uring the lifetime of a cl ump, and v c \ <S v g . 

lLamb and Mi ller (2003) proposed the third version to 
explain that the frequency separation is close to z'spin 
in some stars but close to v sp i n /2 in others. The up- 
per QPO is likewise close to the Keplerian frequency 
vk at r sp . They showed that magnetic and radiation 
fields rotating with the star will preferentially excite ver- 
tical motions in the disk at the 'spin-resonance' radius 
r sr where z/k — z/ s is equal to the vertical epicyclic fre- 
quency. There are two cases in this model. Case 1 sup- 
poses that the flow at r sr is relatively smooth, then the 
vertical motions excited at r sr modulate the X-ray flux 
at v\ w V2 — Vs- Case 1 is fully co mpatible with the 
second version ([Lamb and M iller 200 ll). Case 2 assumes 
that the flow at r sr is highly clumped. In this case, the 
vertical motions excited at r sr modulate the X-ray flux 
at v\ ~ V2 — ^s/2. 

Figure Q] (top panel) displays the fitting result for 4U 
1636-53. We set v\ s s v-i — v s /2 as 4U 1636-5 3 belongs 
to the second case in lLamb and Miller! (]2003l) . The ra- 
tio v c \/vg is represe nted by a free paramete r, the torque 
coefficient cn (see lLamb a nd Miller 120011 for detail). 
Our fitting result is M = 1.545 M©, v s = 652 ± 5 Hz, 
ctv = 0.00274. The fitting does not give a spin frequency 
close to 581 Hz. Moreover, our measured Av ranges from 
220 to 340 Hz, which could be smaller or larger than v s /2 
(2 90.5 Hz). In fact, su ch behavior of Av is already shown 
in Uonker et all (|2002l) . Though successful in explaining 
a changing Av close to half of the spin frequency, the 
latest version predicts that Av is always smaller (larger) 
than z/ s /2 for a prograde (retrograde) flow. When vi is 
below about 800 Hz, Av is larger than 290.5 Hz; other- 
wise, it is smaller than 290.5 Hz. Finally one can notice 
that the fitting result by fixing v s = 581 Hz shows larger 
deviations from the data points. 

The fitting result for Sco X-l is shown in Figure Q] 
(bottom panel), giving M = 1.80 M Q , v s = 329.5 Hz 
(or v s — 659 Hz) and c/v = 0.00216, depending on which 
case we choose. Comp arin g with the fitting resu lt in 
lLamb and Milled (|2001l) and lLamb and Miller! (|200lD , we 
get a larger M and a slightly smaller v s . The reduced % 2 
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Fig. 1. — Fitting results to the sonic-point and spin-resonance 
model for 4U 1636-53 (top) and Sco X-l (bottom). The data points 
with error bars are the observed frequency relations and the curves 
represent the model predictions. In the top panel, the solid curve 
represents the best fitting result, the dashed one shows the fitting 
result by fixing u B = 581 Hz. 

is slightly larger, partly due to our more precise results 
with smaller error bars. 

3.2. Comparisons with orbital resonance models 

iKluzniak and Abramowiczl (|2001j ). lAbramowicz et all 
(2003) introduced several kHz QPO models based on the 
idea of the resonances between the radial and vertical 
frequencies in orbital motion. 

3.2.1. The 2 : 3 parametric resonance model 

A parametric resonance instability occurs near uj r — 
2ojg/n for n = 1,2,3, ••• in an oscillator that obeys a 
Mathieu-type equation of motion, 



; + u)j[l + hcas(oj r t)]56 = 



(3) 



where 68 is the small deviation of elevation 9, the dot 
denotes the time derivative, h is a known constant. 
oj r = 2-KV r and ug — 2~KVg, where v r and vg are 
the radial and vertical epicyclic frequency, respectively 
(|Abramowicz et alll2003D . 



The model predicts v r : vg — 2 : n. When n has 
the smallest possible value, the strongest resonance is 
excited. Since v r < u$, the smallest possible value 
for resonance is n — 3 , meaning that v r : vq = 
2 : 3. Simply supposing v\ = v r and v-i = vg 
(jKluzniak and Abramowiczll2002t ). one can infer the 2 : 3 
ratio of the twin kHz QPO peak frequencies. The ex- 
citation of the r esonance has been stud ied with numeri- 
cal simul ations (lAbramowicz et alll2003l ) and an analytic 
method (|Rebuscoll2004l) . 
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Fig. 2. — Fitting results to the parametric resonance model 
for 4U 1636-53 (top) and Sco X-l(bottom). The dashed curve 
represents the model with v\ : v% = 2 : 3. The dotted curves 
denote the simple assumption (u2 = vg, v\ = u r ) with different M 
an d ?' . The solid curve show s a linear frequency relation according 
to Abramowicz et al (2005). 



The behavior of the frequency relation in the paramet- 
ric resonance model is shown in Figure [2] in comparison 
to our measured data. At first, one can notice that the 
linear 2 : 3 frequency relation (dashed) is not in agree- 
ment with the observations. Then we also notice that 
the model with v 2 = vg and V\ — v r deviates from data 
much more (dotted); it is even unable to follow the ba- 
sic downward-bending track of the data points. Actually 
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the assumption v\ — v r is inappropriate for a NS. Un- 
der the condition of M € [1.4,2.4] M© and j g [0,0.3], 
theoretical v r has a maximum va lue about 635 Hz when 
M = 1.4 Mq, j = 0.3 (see e.g. IStella and Vietril 119991 . 
the equations of orbital frequencies). However, the mea- 
sured v\ reaches as large as 800 ~ 1000 Hz. As shown 
in the figure, the rightmost dotted curve represents the 
predicted frequency relation with M = 1.4 M@, j = 0.3 
and the leftmost one with M = 2.4 Mq, j = 0. The pre- 
dicted curves with other values of M and j lie between 
them. Finally a linear freque ncy relation, i.e. v 2 = v g 
and v 2 = kvi+b, is proposed in lAbramowicz et all ( 2005] ). 
We use it to fit the data and get k = 0.840, b = 395 Hz for 
4U 1636-53, and k = 0.805, b = 422 Hz for Sco X-l, re- 
spectively. Then we get the ratio of v\ to v<i for these two 
sources to be 0.725 and 0.683, respectively, close to but 
higher than 2 : 3. A dopting a function as v 2 = k*v^+b, 
iBelloni et all (|2005t ) also concluded the similar ratios. We 
will discuss the ratios in Section |4] Here we just show 
that the linear fit is naturally disfavored by the data 
points with apparent non-linearity. 

3.2.2. The forced 1 : 2 and 1 : 3 resonance model 

In the numeric al simulations of oscillat ions of a per- 
fect fluid torus (jAbramowicz et a 1 120031) . there is an 
evident resonant forcing of vertical oscillations. The 
forcing is caused by the radial oscillations through a 
pressure coupling. Th is result supports another possi- 
ble resonance model (lAbramowicz and K luzniak 2001: 
lAbramowicz et al 2004). In the model, the resonances 
occur in a forced non-linear oscillator, 

59+LdgS8+[non linear terms in 59} — h(r) cos(uj r t), cug = n 

(4) 

here again u r — 2nv r and luq = 2irvg. 

This model predicts that one of the combination fre- 
quencies, i.e. v- — vg — v r and v + — vg + v ri has a 2 : 3 
ratio to the vertical frequency. For n = 2, the forced 
epicyclic resonance v r : vg = 1 : 2, 



V\ = Vg, 



v 2 = v+, 



(5) 



and for n = 3, the forced epicyclic resonance v r : v$ 
1:3, 



V\ =V-, 



V-2 = Vg. 



(6) 



We fit the observed QPO frequency relations with the 
forced resonance models in Figure [3] The forced 1 : 2 
model predicts that v 2 climbs up to a maximum value at 
v\ w 950 Hz, then decreases rapidly. However, our anal- 
ysis results for the two LMXBs do not show the trend 
that Vi should decrease as v\ increases. Regarding the 
forced 1 : 3 model, it cannot adequately describe the data 
points, especially at the high and low frequencies. The 
vi predicted is higher than that observed at low frequen- 
cies, but lower than that observed at high frequencies. It 
means that the observed Av does not decrease so sharply 
as the model predicts. In addition, these two forced mod- 
els give very large reduced \ 2 . 

3.3. Comparisons with precession models 

This section investigates the predictive ability of three 
precession models, namely, relativistic precession model, 
vertical precession model and total precession model. In 
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Fig. 3. — The fitting results to the forced resonance models 
for 4U 1636-53 (top) and Sco X-l(bottom). The solid curve and 
the dashed one show the forced 1 : 2 and 1 : 3 resonance model, 
respectively. 

these precession models, QPOs can be excited by vari- 
ous resonances with the precession frequencies and or- 
bital frequencies under certain conditions, such as in- 
homo geneities orbiting the inner disk boundary (jStellal 
l200l . 

For the well-known relativistic precession model 
(jStella and "Vietril I1999T ) . the upper QPO v 2 is assumed 
to be the azimuthal frequency the lower QPO V\ 
is expressed as the relativistic periastron precession fre- 
quency, 

v x = v^- v r , (7) 
v 2 = v^. (8) 

In the vertical precession model ([Bursal I2005t ). Vi is 
same as that in Eq. v 2 is hy pothesized as vg . 

In the total precession model (Stuchh'k et al 200^), v\ 
is the total precession frequency, and v 2 is ascribed to the 
Keplerian frequency vk (or the vertical frequency vg), 

V\ = V g - v r (9) 
v 2 = v K or v 2 = vg. (10) 
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Fig. 4. — The fitting results to the precession models for 4U 
1636-53 (top) and Sco X-l (bottom). The solid, dotted and dashed 
curves represent the predictions of the relativistic, vertical and total 
precession models, respectively. The three curves almost overlap, 
especially the relativistic and vertical models. 

Figure 2] illustrates the comparison of the precession 
models with the observed data. Notice that the models 
almost overlap and have the same deviation as the forced 
1 : 3 resonance model: v 2 is predicted too low at high 
frequencies and too high at low frequencies. Hence the 
deviation from the observations increases significantly at 
high and low frequencies. In addition, the relativistic and 
vertical precession models give large M and j. Though 
the total precession model give smaller M and j, the fit 
has largest x 2 among the precession models. 

3.4. Comparisons with disk oscillation models 

The resonances between specific modes in an accretion 
disk are also studied for exciting the observed kHz QPOs. 
In this section, we investigate two models of this kind. 

3.4.1. The deformed-disk oscillation model 

iKatol (|2001l ) brought forward the deformed-disk reso- 
nance model. kHz QPOs are excited by a horizontal res- 
onance in a deformed (warped or eccentric) disk under 
inviscid and adiabatic perturbations. The perturbations 



vary as exp[i(u;i — m<p)], where u is the frequency of the 
perturbations and m (= 0, 1, 2, • • •) denotes the number 
of arms in the azimuthal direction. Various modes o f per- 
turbations are con sidered in a series of their works (|Katol 
I2003L 120051 12009( 1. In the model, QPOs are inertial- 
acoustic oscillations (p-mode) and the gravity oscillations 
(g-mode), or their combination. The twin kHz QPOs are, 



v\ = 2(v K - u r ) 
v 2 = 2vk - v r . 



(11) 
(12) 
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Fig. 5. — The fitting results to the deformed-disk resonance 
model for 4U 1636-53 (top) and Sco X-l(bottom). The solid curve 
represents the best fitting result with M f» 2.4 Mq, j = 0. The 
predicted M and j are almost identical in the two LMXBs. 

Figure [S] exhibits the best fittings for 4U 1636-53 and 
Sco X-l. The fitting re sults (M w 2.4 M , j = 0) are 
consistent with that in Kato (2007). The model gives 
a large NS mass. At the same time, the fittings to the 
two different NS systems give the same set of M and j. 
However the model describes the measured data points 
relatively better than most of other models, though it 
predicts v 2 slightly larger than that observed at high fre- 
quencies. 
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3.4.2. The '-lr, -2v' resonance model 

Unlike the deformed disk oscillation mod el, the per- 
turbations in the '-lr, - 2v' resonance model (jTorok et all 
l2007t lBakala~eFall[2008l) are not stressed in the azimuthal 
direction. In this model, the kHz QPOs are excited by 
the resonance between the radial m = 1 and the vertical 
2 modes. The excited QPOs are supposed to be, 

= v K - i/ r , (13) 

2vk-v 6 . (14) 
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Fig. 6. — The fitting results to the '-lr, -2v' resonance model 
for 4U 1636-53 (top) and Sco X-l(bottom). 



The fitting results to the model are displayed in Figure 
[6l Like the precession and forced 1 : 3 resonance models, 
the model predictions cannot reproduce the observations, 
especially at low and high frequencies. The fittings have 
very large \ 2 an d give M = 2.4 Mq reaching the upper- 
limit in the fitting. 

3.5. Comparison with the higher-order nonlinearity 

model 



Mukhop adhvavl ()2009[ ) treated accreting systems as 
damped harmonic oscillators. These oscillators exhibit 



epicyclic oscillations with higher-order nonlinear reso- 
nance. The resonance is expected to be driven by the 
coupling between the strong disturbance from a NS and 
the weaker one from the flow. In the model, the lower 
and upper kHz QPOs are proposed to be, 



V\ = Vg 



u 2 



n 

V r + ~V S 



(15) 
(16) 



In the disk around a NS, n = 1 corresponds to a 
nonlinear coupling, resulting in Av — v 2 — V\ ~ v s /2\ 
whereas n — 2 corresponds to a linear coupling, result- 
ing in Av i> s . However, n = 3 may correspond to the 
higher-order coupling which is expected to be too weak 
to produce any observable effects. For n = 1 and n = 2, 
the model divides NSs into fast and slow rotators. 

To compute the QPO frequencies, the spin parameter 
j should be determined in the following way. If a NS is 
considered to be spherical in shape with equatorial radius 
R, spin i/ s , mass M, radius of gyration Rq, then the 
moment of inertia and the spin parameter are computed 
by, 

I = MR%, j = JShz, (17) 



GM 2 /c' 



where Q s — 2'kv s . It is known that for a solid sphere 
R 2 G = 2i? 2 /5 and for a hollow sphere R 2 G = 2R 2 /3. How- 
ever a very fast rotating NS is expected to be ellipsoidal 
and not completely solid. Therefore 0.35 < (Rg/R) 2 < 
0.5 is chosen in the model. 

The model parameters, i.e. M, 24, R and (Rg/R) 2 , 
can be obtained by the fitting. For 4U 1636-53, 
Mukhopadhyayl (|2009[ ) treated it as a fast rotator with 
v s = 581 Hz and fit the frequency relation under n = 1. 
In his work, only six data points were collected in the di- 
agram of Ai/ versus z/i . Then he excluded the data points 
at low frequencies, corresponding to the ones with large 
deviations from the model. Finally he gave a low NS 
mass about 1.2 ~ 1.4 Mq. For Sco X-l, the fitting was 
done both under n = 1 and n — 2 and he argued that 
Sco X-l is a slow rotator with n = 2 and v s about 300 
Hz. 

Our fitting results with different n are shown in Figure 
[71 Here we fit all the data points in the two NS systems 
without any exclusion. For 4U 1636-53, the best fitting 
value is M = 1.40 M©, v s = 489 Hz, R = 23.1 km, 
[Rg/R) 2 = 0.38 under n = 1 and M = 1.40, v s = 307 
Hz, R = 26.5 km, (Rg/R) 2 = 0.46 under n = 2. We 
find that under n = 1 and n — 2 the model gives the 
curves almost superposed (solid and dashed in the fig- 
ure). Moreover, v 2 predicted by the model is too high at 
low frequencies. By setting v s = 581 Hz and n = 1, we 
obtain the results with much larger deviations (dotted 
curve). For Sco X-l, the model curves under n = 1 and 
n — 2 almost overlap, resulting in 467 Hz and 294 Hz spin 
frequencies. The best fitting value of other parameters 
is M = 1.40 Mq, R = 23.0 km, (Rg/R) 2 = 0.40 under 
n = 1 and M = 1.40, R = 27.3 km, (Rg/R) 2 = 0.45 
under n = 2. For n = 2, our results are consistent with 
that in lMukhopadhvavl ((20091) . 

3.6. Comparison with the tidal disruption model 
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Fig. 7. — The fitting results to the higher-order nonlinearity 
model for 4U 1636-53 (top) and Sco X-l(bottom). Top panel: the 
solid curve represents the best fitting result under n = 1. The 
dashed curve, which almost overlaps the solid one, is the fitting 
under n = 2. The dotted curve and dot-dashed one denote the 
fittings by setting u B = 581 Hz under n = 1 and n = 2, respectively. 
Bottom panel: the solid and dashed curves indicate the best fitting 
results under n = 1 and n = 2. 

Tidal disruption of the orbits of low-mass satellites 
around a Schwarzschild black hole has recently been 
studied by ICadez et all (j2008l ) . In the clumps of mate- 
rial orbiting such a black hole, a spherical blob can be 
squeezed and stretched by tidal forces into a ring-like 
shape along the orbit, and thus making radial oscillations 
(|Germana et all [20091) . With simulations of such accre- 
tion processes, they generated simulated light curves and 
fit the power spectra of the light curves. Both twin kHz 
QPOs are found and the peak frequencies are supposed 
to be, 

V\ = UK, (18) 

V2 = vk + v r . (19) 

Figure [8] shows the best fittings to this model in the 
two NS systems. As we see, the model describes well the 
main parts of the frequency relations, particularly at low 
frequencies {y\ < 800 Hz). At high frequencies, however, 



Fig. 8. — The fitting results to the tidal disruption model for 
4U 1636-53 (top) and Sco X-l (bottom). The curves exhibits the 
disagreement with the data points when v\ > 800 Hz. 

the model predicts maximum value of V2 and then a sharp 
decrease. It is not supported by observations. Another 
incompatibility is a high NS mass predicted, which is up 
to 2.4 Mq in this model. 

3.7. The Ray leigh- Taylor gravity wave model 

lOsherovich and Titarchukl (|1999l) and iTitarchukl 
(|2003l) described QPOs by the Rayleigh-Taylor insta- 
bility associated with Rossby waves and rotational 
splitting. Twin kHz QPOs are explained as oscillations 
of large scale inhomogeneities (hot blobs) thrown into 
the NS's magnetosphere. Participating in the radial 
oscillations with the Keplerian frequency z^k, such 
blobs are also simultaneously under the influence of the 
Coriolis force. For such mode of oscillations, v-i and 
hold an upper hybrid frequency relation: v\ — v\ = 4r^, 
where v m is the rotational frequency of the magneto- 
sphere near the equatorial plane. If the magnetosphere 
corotates with the NS (solid-body rotation), then the 
spin rotation of the NS would be determined. For the 
first order approximation, v m = v s = const. Within 
the second-order approximation, the slow variation of 



function of vk reveals the structure of the 

magnetospheric differential rotation. Hence, in the 
model, 

ui = uk (20) 



V2 = (v K + 4<) 



2 sl/2 



(21) 



and within the dipole-quadrupole-octupole approxima- 
tion of the magnetic field, the rotation frequency of the 
magnetosphere is, 

Mkk) = Co + Ci4 /3 + C 2 »T + Cs4, (22) 
where Co = v s , C 2 — 2^/C\C^. 
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Fig. 9. — The fitting results to the Rayleigh-Taylor gravity wave 
model for 4U 1636-53 (top) and Sco X-l (bottom). The solid curves 
are the best fitting results while the dash curve in the top panel is 
the fitting result by setting u B = 581 Hz. 

Our fitting treats M, Co, Ci, C3 as free parameters. 
The fitting results are shown in Figure[9] For 4U 1636-53, 
the best fitting (solid curve) returns M = 1.78 M©, Co = 
371, Ci = -0.050, C 3 = -7.7. The spin frequency v s = 
371 Hz, not close to 581 Hz. When we fix v s — 581 Hz, a 
concave curve (dashed) is obtained, whereas the track of 
data points bends downward with a convex shape. For 



Sco X-l, the best fitting values are M = 1.66 M©, Co = 
350, d = -0.046, C 3 = -10.5, respectively. Notice that 
the N S spin frequency of 350 Hz is con sistent with 345 
Hz in lOsherovich and Titarchuk! (|1999|) . Moreover, the 
resulting N S mass 1.66 is reasonab le based on the 
EOS of NS (lLattimer and PrakashlfeOOTj ). 

3.8. Comparisons with the models including the effect of 
magnetohydrodynamics 

The last two models above have included the effect of 
magnetohydrodynamics (MHD) around a rotating NS. 
In a LMXB containing a magnetized NS, the material 
in the accretion disk first rotates in a Keplerian motion, 
then corotates with the magnetosphere as it is trapped 
by the NS magnetic field at the magnetospheric radius, 
and finally flows along the field lines to the polar cap 
of the NS. Some resonant modes may be e xcited by the 
perturbations at the magnetospheric radius ( Zhangl l2004t 
IShi and Lil 12001 . 

3.8.1. The MHD Alfven wave oscillation model 

IZhand (|2004D explained the twin kHz QPOs with the 
MHD Alfven wave oscillations excited by the distortion 
of the NS magnetosphere. The model assumes that the 
infalling MHD material of the Keplerian accretion flow 
distorts the magnetosphere in the regions with enhanced 
mass density gradients, leading to resonant shear Alfven 
waves. In this model, the upper frequency QPO is the 
Keplerian orbital frequency, 

v 2 = vk = 1850AA 3/2 Hz, (23) 

with the parameters X — R/r and A = (m/i?!) 1 / 2 , 
where R 6 = i?/10 6 (cm) and m = M/M Q are the NS 
radius R and mass M in units of 10 6 cm and solar 
masses, respectively. The quantity A 2 is proportional 
to the average mass density of the NS, expressed as, 
(p) = 3M/(4irR 3 ) w 2.4 x 10 14 (g/cm 3 )(^l/0.7) 2 . 

The lower frequency QPO is identified as the Alfven 
oscillation frequency, given as, 



v x = v 2 X 3/i \ll - Vl-X. 



(24) 

Since X is eliminated in the fitting process, we only 
have one parameter A. The comparisons between the 
model predictions and the observations are shown in Fig- 
ure [10] Similar to the precession models, this model also 
predicts v 2 too high at low frequencies and too low at 
high frequencies, leading to the increased deviations from 
the observations at low and high frequencies. The result 
also indicates that Av predicted by this model decreases 
too sharply compared to the observations. The fitting 
for Sco X-l is somewhat better than that for 4U 1636- 
53. Our result of A « 0.7 agrees with that obtained by 
IZhang et all ([2008) , in which the relation of Av versus v 2 
was fitted and the result shows a discrepancy with the 
observations. 

3.8.2. The MHD model 

IShi and L i (2009) presented another explanation for 
kHz QPO signals in LMXBs based on MHD oscillation 
modes in a NS's magnetosphere. Several MHD wave 
modes are derived by solving the dispersion equations. 
They proposed that the coupling of the two resonant 
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Fig. 10. — The fitting results to the MHD Alfven wave model 
for 4U 1636-53 (top) and Sco X-l(bottom). 

MHD modes may lead to the twin kHz QPOs. Finally, 
they presented the following linear frequency relations, 

u 2 = y/l + S 2 ^ + 14), (LSCS) (25) 

^2 = J— g fr* + (SSCS) (26) 

where S 2 = (A 2 -r/ 2 )/(l + ry 2 ) and e 2 = (?? 2 -A 2 )/(l + A 2 ). 
Here A and 77 are two constants linking the Alfven ve- 
locity, acoustic velocity and Keplerian velocity of MHD 
wave in the model. The model divided the twin kHz 
QPOs into two groups with the slope of vs. vi/v s 

relation either larger or smaller than 1.0, i.e., the large 
slope coefficient sources (LSCS) and the small slope co- 
efficient sources (SSCS), respectively. With our fitting, 
we find that both 4U 1636-53 and Sco X-l belong to the 
latter group, since Eq. ([25]) gives bad fitting results with 
reduced \ 2 > 10 5 . 

The fitting results are plotted in Figure [TT] Firstly, 
both for 4U 1636-53 and Sco X-l, the observed relations 
of 1/2 with V\ are not linear; the data points have a track 
to bend downwards. Secondly, for 4U 1636-53, the best 



Fig. 11. — The fitting results to the MHD model for 4U 1636-53 
(top) and Sco X-l(bottom). The solid curves are the best fitting 
results. In the top panel, the dashed curve is the result with v s = 
581 Hz fixed. 

fitting predicts e = 0.65 and v s — 470 Hz. The spin fre- 
quency is not close to 581 Hz. Also, iShi and Lil (|2009fl 
cannot fit well 4U 1636-53 by holding v s = 581 Hz. Their 
result gives e = 0.77 and a large reduced x 2 = 23. In the 
case of Sco X-l, our fitting gives v s = 525 Hz, consider- 
ably larger than that from other models. 

4. DISCUSSION AND CONCLUSION 

We have presented the newly obtained and more ac- 
curate results on the frequency relations of the kHz twin 
QPOs for 4U 1636-53 and Sco X-l. The peak frequen- 
cies of lower and upper QPOs are almost directly pro- 
portional to each other. The data points tend to bend 
downwards between 500 to 1250 Hz in the diagram of 
v\ versus V2 ■ Both of the frequency relations show some 
subtle structure around v\ w 800 Hz. 

Based on the frequency relations, we have systemat- 
ically investigated the predictive ability of all currently 
available models, with which the frequency relation can 
be calculated. Our conclusions are as follows: 

(1) The sonic-point and spin-resonance model seems to 
be only suitable for Z-sources. The model describes well 
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the frequency relation for Sco X-l. However, in the case 
of 4U 1636-53, the model predicts that Av is always less 
than half of the spin frequency for a prograde flow, while 
the observed Av can be larger and smaller than v s /2. 
It indicates that a single rotational direction of accretion 
flow cannot explain the behavior of Av. According to the 
model, such behavior of Av may be produced by a flow 
that retrogrades when it is far away from the NS, but 
then switches to a prograde orbit at some special radius 
when it is closer to the NS. However no evidence is found 
to support such sudden switch of the orbital motion. 

(2) The 2 : 3 parametric resonance model predicts a 
frequency relation bending upwards. Within the limit of 
NS parameters for reasonable NS EOS, it cannot give v\ 
as high as obs ervations. Further more, as claimed i n pre- 
vious papers (|Belloni et al 2005: lZhang~ t al 2006), this 
model leads to the predicted Av increase with the QPO 
frequency. However, the observed downward-bending 
track in the diagram of v\ versus v-i indicates a rough in- 
verse proportion between Av and v\ (or v-i). As regards 
the forced resonance model, the 1 : 2 model predicts a 
decrease of V2 at high frequencies; the 1 : 3 model pre- 
dicts higher (lower) than the observed values at low 
(high) frequencies. In fact, all of the orbital resonance 
models are introduced to explain the observ ed clustering 
of 2 : 3 ratios between v\ and v-i- However. iBelloni et all 
(2005) have demonstrated that a simple random walk 
of the QPO frequencies can reproduce qualitatively the 
observed distributions in freque ncy and frequency ra- 
tio. Later iBoutelier et all (|2009b| ) have pointed out that 
the clustering originates naturally from the sensitivity- 
limited observations and does not support preferred fre- 
quency ratios in NS systems. Our results therefore sug- 
gest that the orbital resonance models should be further 
investigated in order to improve their predictive power 
for the frequency relation. 

(3) All the precession models nearly overlap with each 
other. Their predicted V2 is higher (lower) than that ob- 
served at low (high) frequencies. The deviations from the 
observations increases significantly at high and low fre- 
quencies. In more detail, the relativistic and vertical pre- 
cession models predict a NS mass higher than ~ 2.2 Mq. 
As can be found from Table 3, when we relax the fitting 
limits of M e [1.4, 2.4] and j e [0, 0.3], the extended anal- 
ysis shows that these two models could describe the data 
points better with relatively higher M and j. Essentially, 
the inferred high NS mass may be arisen from the as- 
sumption of the vacuum circumstance around the NS in 
intro ducing the periastron precession term (jZhang et all 
2009). It should be mentioned that considering a small 
eccentricity (< 0.1) which decreases with increasing v$, 
the relativistic precession model would explain the fre- 
quency relation better (Stella and Vietri 1999). In this 
paper, we do not consider the effect of eccentricity on the 
frequency relation because we focus on the NS proper- 
ties. For the total precession model, it gives lower values 
of M and j but larger reduced \ 2 ■ 

(4) The deformed-disk resonance model describes the 
observations relatively better than most of models. The 
fittings give the same M and j for the two LMXBs. It 
suggests that the model may reveal some common prop- 
erties of Atoll and Z sources. Nevertheless the high mass 
predicted and the deviations at high frequencies show 
the model could be modified. For example, the effect 



of magnetic field could be taken into account. For the 
'-lr, -2v' resonance model, it behaves like the procession 
models. Considering that the best fitting results of these 
two models approach the upper limit of M, we also per- 
formed the extended fitting. Both for 4U 1636-53 and 
Sco X-l, the fitting results do not improve significantly. 

(5) The higher-order nonlinearity model classifies NSs 
based on the values of n. After the investigation, one 
can notice that the model predict the nearly identical 
frequency relations under all values of n taken here. 
Thereby, given that we do not know the spin frequency 
of Sco X-l, the classification that Sco X-l is a slow ro- 
tator with n = 2 is not well founded. The superposed 
fitting curves under n = 1 and n — 2 for 4U 1636-53 
also indicate some kinds of ambiguity of the classifica- 
ti on. Apart from t hat, w hen n = 1 is chosen, like that 
in iMu khopadhvav (2009[), the spin frequency predicted 
is not close to 581 Hz. Besides, the model predicts a very 
low NS mass reaching the lower limit in our fitting. Our 
extended analysis shows that the NS masses for the two 
sources are predicted down to 1.0 Mq, quite low for most 
known NS EOS. 

(6) The tidal disruption model can describe the main 
part of the observed frequency relations, though it pre- 
dicts a high NS mass. Maybe it is due to the essential 
difference between NS and black hole systems. This im- 
plies that the kHz QPOs should be greatly affected by the 
surroundings close to the central object. After all, the 
agreement with the observations at low frequencies is re- 
markable, corresponding to the place relatively distant 
to the center. At that place, the clumps of material (or 
particles) are not being exposed to the compact object so 
much. As can be found from Table 3, the fitting results 
of the model can be improved greatly in our extended fit- 
ting, in favor of very large NS masses. The model should 
be modified to better describe the data points at high 
frequencies and to get a more reasonable mass for NS. 

(7) The Rayleigh- Taylor gravity wave model can follow 
the frequency relation in Sco X-l, while for 4U 1636-53 
it cannot predict the 581 Hz spin frequency. This may 
be because the dipole-quadrupole-octupole approxima- 
tion is not sufficiently accurate for the magnetic field. 
Therefore, this model seems promising for explaining the 
origin of kHz QPOs if its description of NS magnetic field 
is more accurate. For the sake of comprehensiveness, it 
should be noted that this model predicts not only the 
high-frequency QPOs, but also the low frequency ones. 
When the low QPOs are also considered, the fittings 
do not always work, unless for one particular frequency 
range one of the low-frequency QPOs is assumed to be a 
harmonic of an unseen one, whereas in the other intervals 
it is the fundamental frequency. 

(8) The MHD Alfven wave oscillation model has the 
same problem as the precession models with increased 
deviations from the observations at high and low frequen- 
cies. It should be noted that the model was put forward 
based on the analogy of the solar coronal atmosphere to 
a NS system. Though the solar coronal atmosphere has 
been studied a lot, the mechanism of Alfven wave os- 
cillations in a NS system remains unclear. The model's 
performance in our fittings indicates that such mecha- 
nism should be investigated further. 

(9) The MHD model predicts a linear frequency rela- 
tion, which is inconsistent with the measured frequency 
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relations. At the same time, the model cannot predict 
reasonable spin frequencies for the two NSs. 
Generally, we also find that: 

(10) These models diverge strongly in their predictions 
of the NS properties. Different models predict spin fre- 
quency from less than 300 Hz to more than 600 Hz. The 
angular momentum is predicted as from to 0.3, cov- 
ering entirely the range of the limit in the fitting. The 
predicted NS mass from different models also covers the 
whole range [1.4,2.4] M . 

(11) The problem of increased deviations at high and 
low frequencies exists in six models: the forced 1 : 3 
model, the three precession models, the '-It, -2v' res- 
onance model and the MHD Alfven wave oscillation 
model. The first five models almost overlap in the plot 
of their fitting results. Since ~ ve (exactly equal if 
v s = 0), these five models have nearly identical expres- 
sions of v\ and v-i- Those five models form the group of 
the largest \ 2 in Table 2. The MHD Alfven oscillation 
model performs slightly better than those five models, 
despite that its \ 2 remains much larger than the other 
remaining models. All of the six models propose that 
the upper frequency QPO is Keplerian, i.e. v-i = v^. 
It infers that for these models, the interpretation is not 
favored by the data. 

(12) Those models including the effects of magnetic 
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field obtain the best fitting results, such as the sonic- 
point beat frequency and the Rayleigh- Taylor gravity 
wave model. At least, they can depict the frequency 
relation for Sco X-l. 

Finally, one should notice the fact that no model gives 
a statistically acceptable % 2 in the fittings. We argue 
that all the models predicting a linear, power-law or any 
other frequency relation are not fully supported by the 
observations, at least for this two sources. 

After the investigation, we comment that since among 
these models we investigated here, three models of 
them (deformed-disk resonance, tidal disruption and the 
Rayleigh- Taylor gravity wave model) have performed rel- 
atively better than other models, we speculate that a 
model which combines these three models together could 
reveal the physical origin of the observed kHz QPO sig- 
nals. It is worth noticing that each of them still has its 
own problems. 
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TABLE 2 

The fitting parameters of all models for 4U 1636-536 and Sco X-l. The errors have been computed by setting 

A X 2 = 1. 







4U 1636-536 






Sco X-l 






Models 


M (M Q ) 


j or spin (Hz) 


Reduced \ 2 


Dof 


M (M ) 


j or spin (Hz) 


Reduced \ 2 


Dof 


SP SR Case 1 


- 


- 


- 


- 


1.80±0.02 


330+^ 


3.0 


8 


SP SR Case 2 


1 545 +0-046 
J J -0.042 


652±6 


7.5 


15 


1.80±0.02 


6591* 


3.0 


8 


fix spin 


1 74fi +0-008 
' -0.009 


581 


18.9 


15 










Para. Res. 


great 


departure 


from 


frequency 


relation 








linear relation 


2.04±0.01 


U. lyztu.ul 


1 1 


14 


2.09±0.01 


0.07±0.01 


71 


7 


17. ........ 1 1 . Q 

rorcea l.o 
Forced 1:2 


, o,r:+U.(W3 

1.010 _ 004 

o nncr +0.006 
2 095 -0.003 


n nnm +0-001U 
0.0001 _ oooi 

n n +0.005 
uu -0.0 


loD 

28 


10 

16 


i.y / ±u.ui 
2.32±0.01 


o.ooo _ 

n n +0.001 



oUO 

54 


y 
9 


Rel. Pre. 


Z.Jia _ 005 


n so + u -° 

U - JU -0.001 


1 

10D 


16 


2 40 +U (J 
Z ' 4U -0.01 


0.250+.0.001 


244 


9 


Ver. Pre. 
Tot. Pre. 


r, 1cn +0.003 

2-160+ oo4 

1 S1/I+0- 004 
1814 -0.003 


0.300 _+™ 
0.0+_°oT 


158 
186 


16 
16 


2.33±0.01 
i 971 +0.001 
- L - t " 1 -0.002 


0.299±0.001 
Oo+oooi 
u - u -0.0 


230 
306 


9 
9 


Deformed-disk 


2 400 +U - U 

Z.4UU _ ooi 


u - uu -0.0 


31 


16 


2 400 +u '" 

Z.4UU _ ooi 


0+u.ooi 

u ' u -0.0 


447 


9 


'-lr, -2v' Res. 


2 400 +U - U 


n „ a +0.002 
a2d8 -0.001 


154 


16 


2 400 +,UI 

^• 4UU -0.002 


n 179 +0.002 
u - 1,z -0.001 


248 


9 


Hi. non-line. n=l 


1.401±0.001 


488.6 1^ 


85 


14 


140 +o.oi 

l w -0.0 


467 ±1 


86 


7 


fix spin 


2 10 + - 01 


581 


211 


15 










Hi. non-line. n=2 


1.401±0.001 


306 7+ 01 


70 


14 


! 40 +o.oi 

L -0.0 


294±1 


70 


7 


Tidal Disruption 


2 400 +U ' U 
-0 .003 


n i«k +0.001 
°- 166 -0.002 


18 


16 


2 400 +U -° 

Z.4UU_ 002 


045+0-001 

u.uio _ .002 


43 


9 


R-T G. Wave 


1 78+0-4V 
L78 -0.38 


371.1 1^ 


10 


14 


1 66 + UM 


350-4 1^ 


4.4 


7 


fix spin 


1 Sfi+0- 38 
1 - 86 -0.40 


581 


38 


15 










Alfven Wave Res. 


A=0.699±0.002 




81 


17 


A=0.658±0.001 




68 


10 


MHD 

fix spin 


t u.u-il _ .014 
c— (1 noo +0-002 
e-0.y28 0.003 


470.3 1 7 ; 7 
581 


10 
21 


16 
17 


e-0 928 +,UUU 


6351^ 


56 


9 



TABLE 3 



The ff 


rriNG para: 


METERS of the 


EXTENDED A> 


IALYSIS 


i for 4U 163 


-6-536 and Sc 


!0 X-l. 








4U 1636-536 






Sco X-l 






Models 


M (M ) 


j or spin (Hz) 


Reduced \ 2 


Dof 


M(M ) j 


or spin (Hz) 


Reduced \ 2 


Dof 


Rel. Pre. 


2.8 


0.49 


131 


16 


3.07 


0.499 


165 


9 


Ver. Pre. 


2.46 


0.5 


141 


16 


2.65 


0.5 


188 


9 


Deformed-disk 


2.48 


0.00 


22 


16 


2.45 


0.0 


395 


9 


'-lr, -2v' Res. 


3.61 


0.500 


129 


16 


3.75 


0.471 


183 


9 


Hi. non-line. n=l 


1.0 


460 


50 


14 


1.0 


438 


61 


7 


Hi. non-line. n=2 


1.0 


283 


34 


14 


1.0 


282 


59 


7 


Tidal Disruption 


3.33 


0.483 


8 


16 


2.70 


0.19 


26 


9 



